// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtFiveFoldTetSubdivisionFilter.C
// ************************************************************************* //

#include <avtFiveFoldTetSubdivisionFilter.h>

#include <cassert>
#include <float.h>

#include <vtkPoints.h>
#include <vtkUnstructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkCellType.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkFloatArray.h>

#include <avtExtents.h>

// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter constructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
// ****************************************************************************

avtFiveFoldTetSubdivisionFilter::avtFiveFoldTetSubdivisionFilter()
{
}


// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter destructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
//  Modifications:
//
// ****************************************************************************

avtFiveFoldTetSubdivisionFilter::~avtFiveFoldTetSubdivisionFilter()
{
}


// ****************************************************************************
//  Method:  avtFiveFoldTetSubdivisionFilter::Create
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
// ****************************************************************************

avtFilter *
avtFiveFoldTetSubdivisionFilter::Create()
{
    return new avtFiveFoldTetSubdivisionFilter();
}


// ****************************************************************************
//  Method:      avtFiveFoldTetSubdivisionFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
// ****************************************************************************

void
avtFiveFoldTetSubdivisionFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const FiveFoldTetSubdivisionAttributes*)a;
}


// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtFiveFoldTetSubdivisionFilter with the given
//      parameters would result in an equivalent avtFiveFoldTetSubdivisionFilter.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
// ****************************************************************************

bool
avtFiveFoldTetSubdivisionFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(FiveFoldTetSubdivisionAttributes*)a);
}


// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter::CreateStitchCells
//
//  Purpose:
//      Indicates that zone information has changed.
//
//  Programmer: Gunther H. Weber
//  Creation:   Mon Jun 20 16:18:56 PDT 2011
//
//
// ****************************************************************************

void
avtFiveFoldTetSubdivisionFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateDataMetaData();
}

// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter::ModifyContract
//
//  Purpose:
//      This method makes any necessay modification to the VisIt contract
//      to, e.g., request that the index variable is also loaded if it.
//      is not already part of the contract.
//
//  Modifications:
//
// ****************************************************************************

avtContract_p
avtFiveFoldTetSubdivisionFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p out_contract = new avtContract(in_contract);

    if (atts.GetAddComponentInformation())
    {
        if( atts.GetIdVar() != "default")
            out_contract->GetDataRequest()->AddSecondaryVariable( atts.GetIdVar().c_str() );
        if( atts.GetValueVar() != "default")
            out_contract->GetDataRequest()->AddSecondaryVariable( atts.GetValueVar().c_str() );
    }

    mainVar = out_contract->GetDataRequest()->GetVariable();
    return out_contract;
}

// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter::PreExecute
//
//  Purpose:
//      Load contour tree before execution
//
//  Programmer: Gunther H. Weber
//  Creation:   Tue Jun 21 14:13:26 PST 2011
//
// ****************************************************************************

void avtFiveFoldTetSubdivisionFilter::PreExecute()
{
    if (atts.GetAddComponentInformation())
    {
        std::ifstream is(atts.GetContourTreeFilename().c_str());
        if (!is.good())
        {
            EXCEPTION1(ImproperUseException, "Need valid contour tree file.");
        }
        contourTree = Branch::unarchive(is);
        //contourTree->print();
        int numEntries = contourTree->maxId() + 1;
        idToBranchMap = new const Branch*[numEntries];
        contourTree->generateIdToBranchMap(idToBranchMap);
        selectedOrHighlighted = new char[numEntries];
        for (int i=0; i<numEntries; ++i) selectedOrHighlighted[i] = 0;
        for (intVector::const_iterator it = atts.GetSelectedIds().begin(); it != atts.GetSelectedIds().end(); ++it)
        {
            if (*it >= 0 && *it < numEntries)
                selectedOrHighlighted[*it] |= 1;
        }
        for (intVector::const_iterator it = atts.GetHighlightedIds().begin(); it != atts.GetHighlightedIds().end(); ++it)
        {
            if (*it >= 0 && *it < numEntries)
                selectedOrHighlighted[*it] |= 2;
        }
    }
}

void avtFiveFoldTetSubdivisionFilter::addBranchIds(vtkIdType vtxIds[4], vtkDataArray *valArray, vtkDataArray *idArray, vtkDataArray *minIdArray, vtkDataArray *maxIdArray, vtkDataArray *sIdArray, vtkDataArray *selectedArray, vtkDataArray *highlightedArray)
{
    int minVtx = 0;
    double minVal = valArray->GetTuple1(vtxIds[0]);
    int maxVtx = 0;
    double maxVal = valArray->GetTuple1(vtxIds[0]);

    double isoval = atts.GetIsovalue();

    for (int i=1; i<4; ++i)
    {
        double val = valArray->GetTuple1(vtxIds[i]);
        if (val < minVal)
        {
            minVtx = i;
            minVal = val;
        }
        if (val > maxVal)
        {
            maxVtx = i;
            maxVal = val;
        }
    }

    Branch::IdxT idMin = idArray->GetTuple1(vtxIds[minVtx]);
    Branch::IdxT idMax = idArray->GetTuple1(vtxIds[maxVtx]);
    minIdArray->InsertNextTuple1(idMin);
    maxIdArray->InsertNextTuple1(idMax);

    //std::cout << "idMin = " << idMin << " idMax = " << idMax << std::endl;
    if (minVal <= isoval && isoval <= maxVal)
    {
        const Branch *b1 = idToBranchMap[idMin];
        assert(b1->id == idMin);
        const Branch *b2 = idToBranchMap[idMax];
        assert(b2->id == idMax);
        Branch::IdxT isosurfBranchId = Branch::branchIdForValue(b1, b2, isoval); 
        //std::cout << "====== ID for surface between " << idMin << " and " << idMax << " for value " << isoval << " is >>>>>>>>>>>>>>>>>" << isosurfBranchId << std::endl;
        sIdArray->InsertNextTuple1(isosurfBranchId);
        selectedArray->InsertNextTuple1((selectedOrHighlighted[isosurfBranchId] & 1) ? 1 : 0);
        highlightedArray->InsertNextTuple1((selectedOrHighlighted[isosurfBranchId] & 2) ? 1 : 0);
    }
    else
    {
        sIdArray->InsertNextTuple1(-1);
        selectedArray->InsertNextTuple1(0);
        highlightedArray->InsertNextTuple1(0);
    }
}


// ****************************************************************************
//  Method: avtFiveFoldTetSubdivisionFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the FiveFoldTetSubdivision filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Mon Jun 20 14:54:18 PST 2011
//
//  Modifications:
//    Eric Brugger, Mon Jul 28 15:28:00 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtFiveFoldTetSubdivisionFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    vtkPoints *points = vtkPoints::New(); 
    int dims[3];
    if (vtkRectilinearGrid* rgrid = dynamic_cast<vtkRectilinearGrid*>(in_ds))
    {
        rgrid->GetDimensions(dims);
        points->SetNumberOfPoints(rgrid->GetNumberOfPoints());

        for (int ptId=0; ptId<rgrid->GetNumberOfPoints(); ++ptId)
        {
            points->SetPoint(ptId, rgrid->GetPoint(ptId));
        }
    }
    else if (vtkStructuredGrid* sgrid = dynamic_cast<vtkStructuredGrid*>(in_ds))
    {
        sgrid->GetDimensions(dims);
        points->ShallowCopy(sgrid->GetPoints()) ;
    }
    else
    {
        EXCEPTION1(ImproperUseException, "Mesh must be rectilinear or structured!");
        points->Delete();
        return in_dr;
    }

#if 0
    if ((in_ds->GetCellData()->GetScalars()) || (in_ds->GetCellData()->GetVectors()))
    {
        EXCEPTION1(ImproperUseException, "Cell centered data is not supported!");
        points->Delete();
        return in_ds;
    }
#endif

    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints(points);
    points->Delete();

    vtkDataArray *valArray = 0;
    vtkDataArray *idArray = 0;
    vtkDataArray *minIdArray = 0;
    vtkDataArray *maxIdArray = 0;
    vtkDataArray *sIdArray = 0;
    vtkDataArray *selectedArray = 0;
    vtkDataArray *highlightedArray = 0;
    if (atts.GetAddComponentInformation())
    {
        if (atts.GetIdVar() != "default")
        {
            idArray = in_ds->GetPointData()->GetScalars(atts.GetIdVar().c_str());
        }
        else
            idArray = in_ds->GetPointData()->GetScalars();

        if (!idArray)
        {
            EXCEPTION1(ImproperUseException, "Need vertex centered ID data!");
        }
        if (atts.GetValueVar() != "default")
        {
            valArray = in_ds->GetPointData()->GetScalars(atts.GetValueVar().c_str());
        }
        else
            valArray = in_ds->GetPointData()->GetScalars();
        if (!idArray)
        {
            EXCEPTION1(ImproperUseException, "Need vertex centered value data!");
        }
 
        minIdArray = vtkFloatArray::New();
        minIdArray->SetName("MinID");
        maxIdArray = vtkFloatArray::New();
        maxIdArray->SetName("MaxID");
        sIdArray = vtkFloatArray::New();
        sIdArray->SetName("SurfaceID");
        selectedArray = vtkFloatArray::New();
        selectedArray->SetName("Selected");
        highlightedArray = vtkFloatArray::New();
        highlightedArray->SetName("Highlighted");
 
        assert(idArray);
        assert(valArray);
    }

    for (int k=0; k<dims[2]-1; ++k)
        for (int j=0; j<dims[1]-1; ++j)
            for (int i=0; i<dims[0]-1; ++i)
            {
                //vtkIdType cellId = (k*(dims[1]-1) + j)*(dims[0]-1) + i;
#define POINT_ID(I,J,K) ((K)*dims[1]+J)*dims[0]+I
                vtkIdType pt0Id = POINT_ID(i  , j  , k  );
                vtkIdType pt1Id = POINT_ID(i+1, j  , k  );
                vtkIdType pt2Id = POINT_ID(i+1, j+1, k  );
                vtkIdType pt3Id = POINT_ID(i  , j+1, k  );
                vtkIdType pt4Id = POINT_ID(i  , j  , k+1);
                vtkIdType pt5Id = POINT_ID(i+1, j  , k+1);
                vtkIdType pt6Id = POINT_ID(i+1, j+1, k+1);
                vtkIdType pt7Id = POINT_ID(i  , j+1, k+1);
#undef POINT_ID

                int vtx0HasSixNeighborsParity = atts.GetOddParityHasSixNeighborhood() ? 0 : 1;

                if (((i+j+k) % 2) == vtx0HasSixNeighborsParity)
                {
                    vtkIdType tet0[4]  = {pt1Id, pt4Id, pt3Id, pt0Id};
                    vtkIdType tet1[4]  = {pt1Id, pt3Id, pt6Id, pt2Id};
                    vtkIdType tet2[4]  = {pt1Id, pt6Id, pt4Id, pt5Id};
                    vtkIdType tet3[4]  = {pt3Id, pt4Id, pt6Id, pt7Id};
                    vtkIdType tet4[4]  = {pt1Id, pt3Id, pt4Id, pt6Id};
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet0);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet1);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet2);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet3);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet4);
                    if (minIdArray && maxIdArray && sIdArray)
                    {
                        addBranchIds(tet0, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet1, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet2, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet3, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet4, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                    }
                }
                else
                {
                    vtkIdType tet0[4]  = {pt0Id, pt2Id, pt5Id, pt1Id};
                    vtkIdType tet1[4]  = {pt0Id, pt7Id, pt2Id, pt3Id};
                    vtkIdType tet2[4]  = {pt0Id, pt5Id, pt7Id, pt4Id};
                    vtkIdType tet3[4]  = {pt2Id, pt7Id, pt5Id, pt6Id};
                    vtkIdType tet4[4]  = {pt0Id, pt7Id, pt5Id, pt2Id};
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet0);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet1);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet2);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet3);
                    ugrid->InsertNextCell(VTK_TETRA, 4, tet4);
                    if (minIdArray && maxIdArray && sIdArray)
                    {
                        addBranchIds(tet0, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet1, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet2, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet3, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                        addBranchIds(tet4, valArray, idArray, minIdArray, maxIdArray, sIdArray, selectedArray, highlightedArray);
                    }
                }
            }

    ugrid->GetPointData()->ShallowCopy(in_ds->GetPointData());

    if (minIdArray && maxIdArray && sIdArray)
    {
        ugrid->GetCellData()->AddArray(minIdArray);
        ugrid->GetCellData()->AddArray(maxIdArray);
        ugrid->GetCellData()->AddArray(sIdArray);
        ugrid->GetCellData()->AddArray(selectedArray);
        ugrid->GetCellData()->AddArray(highlightedArray);
        //ugrid->GetCellData()->SetActiveScalars("SurfaceID");
        minIdArray->Delete();
        maxIdArray->Delete();
        sIdArray->Delete();
        double range[2] = { FLT_MAX, -FLT_MAX };
        GetDataRange(ugrid, range, "SurfaceID", false);
        avtDataAttributes &dataAtts = GetOutput()->GetInfo().GetAttributes();
        TRY {
            dataAtts.GetThisProcsOriginalDataExtents("SurfaceID")->Set(range);
            dataAtts.GetThisProcsActualDataExtents("SurfaceID")->Set(range);
            range[0] = 0;
            range[1] = 1;
            dataAtts.GetThisProcsOriginalDataExtents("Selected")->Set(range);
            dataAtts.GetThisProcsActualDataExtents("Highlighted")->Set(range);
        }
        CATCHALL
        {
        }
        ENDTRY
    }

    avtDataRepresentation *out_dr = new avtDataRepresentation(ugrid,
        in_dr->GetDomain(), in_dr->GetLabel());

    ugrid->Delete();

    return out_dr;
}
